Representing perturbed dynamics in biological network models 



Gautier StoU^'^Q Jacques Rougemont^j^ and Felix Naef^'^-^ 
^ NCCR Molecular Oncology, ch. des Boveresses 155, 1066 Epalinges, Switzerland 
^School of Life Sciences, ISREC, Ecole polytechmque Federale de Lausanne 1015 Lausanne, Switzerland and 
^ Swiss Institute of Biomformatics, Quartier Sorge-Genopode, 1015 Lausanne, Switzerland 

We study the dynamics of gene activities in relatively small size biological networks (up to a few 
tens of nodes), e.g. the activities of cell-cycle proteins during the mitotic cell-cycle progression. 
Using the framework of deterministic discrete dynamical models, we characterize the dynamical 
modifications in response to structural perturbations in the network connectivities. In particular, 
we focus on how perturbations affect the set of fixed points and sizes of the basins of attraction. 
Our approach uses two analytical measures: the basin entropy H and the perturbation size A, a 
quantity that reflects the distance between the set of fixed points of the perturbed network to that 
of the unperturbed network. Applying our approach to the yeast-cell cycle network introduced by 
Li et al. provides a low dimensional and informative fingerprint of network behavior under large 
classes of perturbations. We identify interactions that are crucial for proper network function, and 
also pinpoints functionally redundant network connections. Selected perturbations exemplify the 
breadth of dynamical responses in this cell-cycle model. 
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I. INTRODUCTION 



Recent experimental developments in the fields of ge- 
nomics, e.g. whole genome DNA sequencing or pro- 
teomics, are opening possibilities for systems level studies 
in biology [U El [3j |4] . In particular, the notion that bio- 
logical functions may rely on a large number of intercon- 
nected variables (for example genes) working in concert 
has stimulated general theoretical interest about prop- 
erties of biological networks [5j. Studies of the statis- 
tical properties of large (typically thousands of nodes) 
biological networks have identified a number of func- 
tional building block, termed network motifs, that occur 
more frequently than random These findings sup- 
port the idea that some systems are designed around a 
modular architecture, in which autonomous modules are 
wired together to generate versatile biological functions 
[TJ m [HI [23 • While structural (or topological) prop- 
erties are key for network characterization, functional 
properties are ultimately encoded in dynamical, or time- 
dependent changes in the state variables of the nodes. 
The sizes of systems that can be modeled dynamically 
are typically much smaller (10-100 nodes). One common 
modeling approach, for example for the yeast cell-cycle 
[lOj . is to simulate the nonlinear system of chemical rate 
equations describing the putative biochemical processes. 
Modeling approaches have been applied to a number of 
systems, including the cell-cycle [101 [II]j the lambda- 
phage switch in E. coli [9 . Although these models pro- 
vide a detailed description, this approach suffers from the 
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caveat that most parameters are currently not accessible 
experimentally. In addition, the number of parameters 
is typically about five per reaction, resulting in a pro- 
hibitively large parameter space. This last point makes it 
difficult to grasp the full solution space of the model. Re- 
cent approaches based on sampling the parameter space 
in optimal regions have been developed [19]. At the op- 
posite end of model complexity, dynamical rules based 
on boolean state variables have been useful for studying 
more global dynamical properties of topological classes 
of networks [20l [21]. In addition, boolean models have 
been successfully applied to the yeast cell-cycle [22, 25] 
and the body patterning in drosophila embryos [23, .24j . 

In this study, we develop a systematic approach to de- 
scribe how the dynamical landscape of small (less than 
about 50 nodes) boolean networks is affected by pertur- 
bations in the network connectivity. In particular, we 
consider the basin entropy H, a quantity that consid- 
ers the size distribution of the basins of attraction. We 
complement entropy with a measure of distance between 
the stable fixed points of a perturbed network and those 
in the unperturbed network. This combination gives a 
low-dimensional and compact representation of the pat- 
terns induced by a large number of perturbations. We 
illustrate our methods using the yeast cell-cycle network 
introduced in |22| . and discuss examples of structural 
perturbations producing a range of modified basins of 
attraction. 



II. DEFINITIONS 

Following [22] a network of N nodes can be represented 
by a iV X adjacency matrix A, in which an activating 
link between node i and node j is represented by Aij — 
1 and an inhibiting link by Aij = —1. The possibility 
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of self-inhibitory (or activating links) An = ±1 is not 
excluded. In the Boolean approximation, each node has 
two possible states, so that the global state of all nodes 
can be represented by a vector S, with Si — I when the 
node i is on and 5^ = if the node is ojf. The full phase 
space containing 2^ states is denoted by A. 

A. Boolean dynamics 

A simple dynamical rule that characterizes the tempo- 
ral evolution of the state variable can be defined following 
[35], which is closely related to update rules applied in 
perceptron models. If the network is in the state S{t) at 
time t, the state at the next time-step S{t + 1) is given 
by: 

r 1 if E,-4„-^,W >0 
S,{t+1)={ S^it) if E,A.jS,{t)^0 (1) 
[ if j:^A,S,{t)<0 

For a given network, we apply this rule to every pos- 
sible initial condition in A. This defines orbits (trajecto- 
ries) that must end in a limit cycle (periodic attractor) 
since we are dealing with a dynamical system on a finite 
space. A fixed point is a cycle of length one. 

Accordingly, A can be decomposed into a disjoint union 
of K basins of attraction Bk of size dfe: A = IJ^-^ 

In a biological network, the attractors correspond to 
functional endpoints, and it is important that the states 
in the attractors are consistent with observed data. For 
example, by far the largest endpoint in the cell-cycle net- 
work of Li et al. (see appendix) corresponds to the sta- 
tionary Gl phase in the cycle. Other systems are more 
switch-like, for instance in signal transduction, where a 
cell might change its state from growth to differentiation 
according to an external trigger. To characterize these 
attractors, we introduce the following definitions: 

• We compute the number of attractors K: an attrac- 
tor is a limit cycle or a fixed point. An attractor A 
has a basin of attraction B which is the set of all 
initial conditions whose orbit converges to A. 

• The basin entropy H is defined as follows: let 
Pk ~ 2^^ dk be the probability that an initial state 
belongs to basin B^. Then, the entropy reads 

K 

H:=-Y,Pklog{pk) (2) 

A;=l 

H is maximum {H — log(_ft')) if each state is its 
own basin of size one, and minimum (H = 0) when 
there is one single basin. iJ is a natural measure 
for characterizing basin structures [27]. Because it 
takes into account the relative basin sizes, it is quite 
insensitive to appearance of small and biologically 
irrelevant basins. 



• The perturbation size A measures the distance be- 
tween attractors of a perturbed and a reference net- 
work: from every initial conditions, the Hamming 
distance between the fixed points is computed, and 
the average over all initial conditions is taken. More 
precisely, if FPg(S) is the fixed point of the tra- 
jectory starting at S and generated by the network 
G, then 

^G,G' := ^5]HAM(FPg(S),FPg'(S)) (3) 
s 

where HAM(-, •) is the Hamming distance between 
two boolean states, namely 

HAM(S,T):=1^|S,-T,| (4) 

i 

The value A has the following interpretation: it is 
the average probability (taken over all nodes) that, 
for a random initial condition, the final state of a 
node differs. In this study, the reference network 
G will be the cell-cycle network of Li et al., which 
has one very a large basin of attraction and several 
smaller ones. If some trajectories in the perturbed 
networks G' end in a limit cycle, A is defined as the 
average of the Hamming distance along the cycle. 

B. Network models and perturbations 

Our goal is to assess how network dynamics is af- 
fected by several types of perturbations. We consider 
two classes: one which randomizes the adjacency ma- 
trix while keeping a number of topological characteristics 
from the original network invariant. The second class 
mimics biological perturbations, as would occur for ex- 
ample through mutations in the interaction partners that 
constitute the network links. The two classes are defined 
as follows: 

• Shujfle (class I) : all activating and inhibiting arrows 
are cut in half and re- wired randomly. This ensure 
that the connectivity at each node is conserved. 
As compared to the Li et al. [32] study, we gen- 
erate random networks that are more constrained, 
since the connectivity at each node is forced to re- 
main unchanged after randomization. Such pertur- 
bations are applied in the studies of network motifs 

mm- 

• Remove (class II): the arrows are simply sup- 
pressed. We extend this class of perturbations be- 
yond single link removal. 

III. RESULTS 

We study the yeast cell-cycle network of Li et al. [33] 

(the Yeast cell-cycle network or YCC), in which a 
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boolean model reproducing the different phases of the 
cycle is constructed (see appendix). This model has a 
main fixed point attracting 86% of the intial conditions. 
Biologically this state corresponds to the Gl stationary 
phase of the cell-cycle, as reflected by the activities of 
the respective nodes. Using computer simulations, the 
authors further showed that the cell-cycle dynamics had 
certain robustness properties when challenged with per- 
turbations. In particular, it was shown that in a majority 
of cases, removal of one link or addition of a link at ran- 
dom did not change much the size of the largest basin of 
attraction. Finally, the studied network had unusual tra- 
jectory channeling properties, when compared to random 
networks with equal number of nodes and links. Here we 
extend the characterization of this model by introducing 
a combination of measures to characterize the structure 
of basins of attraction as they are modified by structural 
perturbations. In particular we investigate the conse- 
quences of combined mutations and show that they can 
lead to cancellation effect. 
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A. Study of shuffled networks (Class I 
perturbations) 

This type of perturbation allows to study the dynam- 
ical characteristics of a biological network in comparison 
with random networks belonging to a topological class. 
Figure 1 shows the Number of attractors (K) and the 
Entropy (H) of the YCC and randomly shuffled (Class 
I) versions thereof. 

The location of the reference network in the H — K 
plane respective to the scatter of the perturbed networks 
allows us to asses how typical a network behaves with 
respect to a class. Accordingly, the YCC is atypical, as 
seen by its marginal location in the lower left corner. 
Indeed, this network has lower entropy and fewer basins 
than most networks, consistent with ^T2\. 



B. Study of mutated networks (Class II 
perturbations) 

The previous discussion shows how entropy character- 
izes the system of attractors. However, H contains only 
information about the relative weights of the attractors, 
irrespective of their biological relevance. For example a 
perturbation can decrease the entropy while shifting the 
fixed point away from from that in the unperturbed, bi- 
ologically relevant state. For this reason we introduced 
a second quantity, A (Equation [s]) , a probabilistic mea- 
sure of the change in the fixed point after perturbation. 
Therefore, A reflects the change in the biological rele- 
vance of the basin structure. 

We first repeat Figure 1 for class II perturbations 
which shows that networks with few perturbations clus- 
ter around the wild-type model (Figure 2A), while the 
sread for networks with four perturbations resembles the 



FIG. 1: Entropy vs. number of attractors after class I per- 
turbation (shuffled arrows). The range of possible H values 
is indicated by the dashed gray lines. The open red circle 
represents the reference network, the other points show the 
perturbed networks. 



shuffled models (Figure 1). Turning to the measure of A, 
we find that A-distribution (Figure 2B) is bimodal, show- 
ing two distinct populations of perturbations: (A < 0.2 
and A > 0.2). In the second case, the perturbed model 
does not reproduce the biologically correct cell-cycle pro- 
gression. But if A is small, then the system of attractors 
of the perturbed network is still consistent with the biol- 
ogy and entropy allows to discriminate between networks 
with a larger or smaller main basin of attraction. For this 
reason, the entropy and A are complementary for de- 
scribing the dynamical landscape (Figure 2C). The two 
different modes in the A-histogram are clearly reflected 
on this 2D representation. Noticeably, the A values span 
a broad range for any number of removed arrows, on 
the other hand higher entropies are more frequent for 
larger number (> 2) of removed arrows. Qualitatively, 
the spread of points in the H — A plane conveys a mea- 
sure of robustness. Accordingly, the A measure appears 
more fragile than the entropy property, especially when 
few arrows are removed. 

We now interpret the different locations in the H — A 
plane: 

1. If A is large (A > 0.2), the model does have attractor 
states which coincide with the gene activities of the 
different cell-cycles phases. Such perturbations are 
specially interesting if the number of removed arrows 
is small (dark colors). Such links are then essential for 
the model, as their removal disrupts the cell-cycle very 
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C. Examples of mutations 

In the first example (Table II), the dynamics has a 
large main basin of attraction like in the unperturbed 
model (Table I). However, the fixed point is significantly 
different fi:om wild-type as the system is blocked in the a 
state of the M-phase and cannot finish properly the cell- 
cycle (see Appendix for the recapitulation of the wild- 
type mode from [H]). 



Cln3 MUF SBF Clnl,2 Cdhl SwiS Cdc20,14 Clb5,6 Sicl Clbl,2 Mcml % 
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0.8613 
0.0737 
0.0532 
0.0043 
0.0034 
0.0034 
0.0004 



TABLE I: Basins of attraction with their respective probabil- 
ities in (%) for the original YCC network. The largest basin 
ends at the Gl stationary state. Entropy H — 0.543, Number 
of attractors K = 7. 

In the second example (Table III), the dynamics has 
the same main fixed point as the wild-type, but with a 
smaller basin of attraction, while the second biggest has 
grown. Therefore the removed connection SBF Clnl,2 
contributes to the ability of the main fixed point to funnel 
trajectories. 

The third example (Table IV) is a model with four re- 
moved arrows which has the same main fixed point with a 
slightly higher probability. Also, the second largest fixed 
point is same as in the wild-type model. This indicates 
that the effect of some mutations can be canceled by fur- 
ther mutations. While such cases exist, we found that 
networks with several removed links that preserving the 
unperturbed cell-cycle behavior are rare. 



FIG. 2: Entropy, number of attractors and A after class II 
perturbation (removed arrows). Colors represent different 
number of removed arrows: black for one removed arrows, 
red for 2, green for 3, turquoise for 4 and yellow for more 
than 4. A: same figure as for the class I perturbation, the 
range of possible H values is indicated by the dashed gray 
lines, the open blue circles represent the reference network. 
B: Distribution of A. C; A vs. H plot, the dashed gray line 
represents the entropy of the reference network. 



efficiently. 

2. If A is small and the entropy increases, the proba- 
bility that the dynamics ends in the reference attractor 
decreases demonstrating that the removed arrows con- 
tributed to the channeling properties of the system. 

3. If A is small and the entropy decreases, the main 
attractor of the perturbed network has a stronger 
attraction property. Some of these networks could be 
considered as alternative cell-cycle models. 

We illustrate these three regimes by examples: 



IV. CONCLUSION 

We have proposed a systematic approach for studying 
the dynamical attractor landscape of biological networks, 
and their response to structural perturbations. In partic- 
ular, we introduced a low dimensional representation of 
the system of attractors, the entropy, and a probabilistic 
measure in the perturbation size A. This enabled us to 
study the global characteristics of network perturbation 
in a compact and visually effective form. In a biological 
context, this can provide hints to elucidate the dynamical 
role of specific network links. Alternatively, the function 



Cln3 MBF SBF Clnl,2 Cdhl SwiS Cdc20,14 Clb5,6 Sicl Clbl,2 Mcml 



0.880 
0.054 
0.027 
0.015 
0.010 
0.004 
0.003 
0.003 
0.000 



TABLE II: Basins of attraction with their respective prob- 
abilities, when (Cdc20,Cdcl4) Clbl,2 and Sicl ^ Clbl,2 
are removed. Entropy — 0.549, Number of attractors = 9, A 
= 0.41. 
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Cln3 MBF SBF Clnl,2 Cdhl SwiS Cdc20,14 Clb5,6 Sid Clbl,2 Mc 



% 

0.6669 
0.1762 
0.0654 
0.0532 
0.0180 
0.0043 
0.0043 
0.0034 
0.0034 
0.0034 
0.0004 
0.0004 
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APPENDIX A: THE YEAST CELL-CYCLE 
NETWORK OF LI ET AL. 



TABLE III: Basins of attraction with their respective proba- 
bihties, when SBF — » Clnl,2 is removed. Entropy H — 1.096, 
Number of attractors K = 12, A ^ 0.05. 



The following two tables are recapitulated from 



Cln3 MBF SBF 01iil,2 Cdhl SwiS Cdc20,14 Clb5,6 Sid Clbl,2 Mcml % 



10— >2 10— '3 10^5 



10— >6 10^9 



.'lO 7^8 7-^10 8^5 8- 
1^1 4^4 6^6 7^7 11- 



0.8793 
0.0507 
0.0356 
0.0268 
0.0034 
0.0034 
0.0004 



TABLE IV: Basins of attraction with their respective prob- 
abilities, when (Cdc20,Cdcl4) Clbl,2, Clbl,2 Mcml, 
Clbl,2 — > Cdhl and Clbl,2 — > Swi5 are removed. Entropy 
H = 0.523, Nmnber of attractors K = 7, A = 0.025. 



of new and yet unobserved links can be predicted as in 
[26j . and imperfect starting models can be improved. 

We applied this method to a model of the yeast cell- 
cycle by Li et al. Using the measures introduced, we 
have generalized the dynamical characterization of the 
model using a broad range of perturbations. This has 
enabled us to emphasize the breadth of dynamical be- 
havior (Figure 2) induced by only few mutated hnks. In- 
terestingly, we observed (Figure 2C) that the structure 
of the system of attractors {H) behaves quite robustly 
compared to the modification in the final states (A), es- 
pecially when the number of removed links is small (< 3). 
We illustrated through examples the consequences of re- 
moving individual or groups of links. Interestingly it was 
possible to remove up to four links while not affecting 
the basin structure significantly. Tracking the dynamical 
changes in the activity levels of proteins in a network is 
a very high-dimensional problem. It therefore important 
to be have few informative variables which allow one to 
efficiently assess a large number of perturbed models at 
once. We believe that basin entropy and distance to a 
reference attractor are well suited for this purpose. 



TABLE V: Adjacency matrix of the Yeast cell-cycle network. 
The numbers refer to the ordering of the nodes as used in 
Tables I-IV,VI. -I- (respectively — ) represent activating (re- 
spectively repressing) links. 



t Cln3 MBF SBF Clnl,2 Cdlil 



C20.14 Clb5,6 Sid Clbl,2 Mcml Pha 



START 
Gl 
Gl 
Gl 
S 

G2 



Gl 
Gl* 



TABLE VI: This table represents the discrete time evolution 
of the boolean states of the YCC network as it traverses the 
different cell-cycle phases. Cdc20.14 has been abbreviated 
C20,14; Gl* indicates the stationary Gl phase. 
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